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Abstract. We have investigated structure and properties of small metal clusters 
using all-electron ab initio theoretical methods based on the Hartree-Fock 
approximation and density functional theory, perturbation theory and compared 
results of our calculations with the available experimental data and the results of 
other theoretical works. We have systematically calculated the optimized geometries 
of neutral and singly charged sodium clusters having up to 20 atoms, their multipole 
moments (dipole and quadrupole), static polarizabilities, binding energies per atom, 
ionization potentials and frequencies of normal vibration modes. Our calculations 
demonstrate the great role of many-electron correlations in the formation of electronic 
and ionic structure of small metal clusters and form a good basis for further detailed 
study of their dynamic properties, as well as structure and properties of other atomic 
cluster systems. 



1. Introduction 

Atomic clusters and small nanoparticles have been recognized as new physical objects 
with their own properties relatively recently. This became clear after such experimental 
successes as the discovery of electron shell structure in metal clusters |J, observation 
of plasmon resonances in metal clusters [@-@] and fullerenes |5|,||, formation of singly 
and doubly charged negative cluster ions @ and many more others. The novelty of 
cluster physics is also greatly connected with the fact that cluster properties explain the 
transition from single atoms or molecules to solid state. Comprehensive survey of the 
field can be found in review papers and books, see e.g. [j3Hl4j]. 

There are many different types of clusters, such as metallic clusters, fullerenes, 
molecular clusters, semiconductor clusters, organic clusters, quantum dots, positively 
and negatively charged clusters, which all have their own features and properties. In 
this paper we focus on the detailed systematic study of the structure and properties 
of small metal clusters and in particular sodium clusters using ab initio all-electron 
many-body theory methods. 
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So far, for sodium clusters, systematic calculations of cluster properties on the same 
level of theory as in our present work (i.e. all electron ah initio) have been performed 



only for clusters with N < 10 [13, 15H19||, where N is a number of atoms in a cluster. 
In our work we extend this limit up to iV < 20. Note that most of the cited papers 
are focused on the investigation of neutral cluster properties rather than ions. In our 
present work we perform systematic comparative analysis of properties of neutral and 
singly charged sodium clusters in the specified size range. 

During the last decade, there were performed numerous experimental and 
theoretical investigations of the properties of small metal clusters as well as the processes 
with their involvement. Here we are not able to review even all essential results 
obtained in the field and only refer to those, which are related the most closely to 
the subject of our paper. In []]], it was experimentally proved that metal clusters have 
the shell electronic structure and the magic cluster numbers have been determined by 
observation of the sodium cluster abundances in mass spectra. Experimental study 
of electronic structure and properties of small metal clusters have been performed 
20|.|2T| (for review also see P,[l0|,p!l|,[l3|,p^|). In [pC(] , there have been measured the 



in 



ionisation potentials for a sequence of small neutral and positively charged sodium metal 
clusters, which independently proved their shell structure. The dipole polarizabilities 
of sodium clusters have been experimentally determined in [[H]]. Dissociation energies 
of neutral and positively charged small sodium and potassium metal clusters have been 
measured in p2|-p4|j. Dynamical properties of clusters have been studied by means of 
photon, electron and ion scattering. These methods are the traditional tools for probing 
properties and internal structure of various physical objects. Using these methods, 
for example, plasmon excitations in metal clusters [@, and fullerenes have been 
observed (for review also see P,pH|]). 

Metal clusters have also been studied theoretically. Structural properties of small 
metal clusters have been widely investigated using quantum chemistry methods. Here 
we refer to the papers Hl5|-[T8|, |26l-f29|[ , in which optimized geometries, binding energies, 
ionization potentials, electron structure and electron transport properties of small 
lithium and sodium clusters have been calculated. In these papers the systematic 
analysis of the cluster properties has been limited by cluster sizes N < 10. In the 
present paper we extend this limit up to iV < 20 and perform systematic analysis of 
various cluster characteristics both for neutral clusters and singly charged cluster ions. 

In a last few years, a number of papers have been devoted to the calculation of 



dipole static polarizabilities of neutral sodium and lithium clusters [19,30-33]. Note that 



most of these studies have been performed within the cluster size range N < 20. The 
results of different theoretical approaches have been compared with the experimental 
data from pi]] . However, only in |B||, calculations of the cluster geometries and 
polarizabilities have been performed on the same level of theory (i.e. all electron ab 
initio) as in our work and were limited by N < 8. 

Alternatively, the jellium model for metal clusters was suggested. This model 
explains well enough the shell structure of metal clusters and their essential dynamic 
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properties, such as plasmon excitations. Initially, jellium calculations for metal clusters 
were based on the density functional formalism with the use of pseudopotentials for 



the description of electron relaxation effects and lattice structure [33]. Fully self- 
consistent calculations for spherical jellium metal clusters have been performed within 
the framework of the spin-density-functional method [55] and the Kohn-Sham formalism 
for the self-consistent determination of electron wave functions |j3T5|, |57[]. The Hartree- 
Fock scheme for the self-consistent determination of the electron wave functions of 
spherical jellium metal clusters was also introduced later in [[3^, |39|. This approach 
was generalized for axially deformed cluster systems in flO| . Dynamical jellium model 
for metal clusters, which treats simultaneously collective vibrational modes (volume 
vibrations, i.e. breathing, plus shape vibrations) of the ionic jellium background in a 
cluster, the quantized electron motion and interaction between the electronic and ionic 
subsystems was developed in [(41~] . |42|1 . 

The jellium model provides a very useful basis for studying various collision 
processes, such as photabsorption [|43|] , photoionization || |44j, |45[] , elastic [0,0 anc ^ 



inelastic scattering p7|-[50|[, electron attachment [^T],^2j, photon emission |53],|54j and 
others, involving metal clusters,. On the basis of the jellium model one can develop 
ab initio many-body theories, such as the random phase approximation with exchange 
or the Dyson equation method and effectively solve many-electron correlation problem 
even for relatatively large cluster systems containing up to 100 atoms or even more. 
Review of these methods in their application to the electron scattering of metal clusters 
one can find in p5|. As elucidated in the papers cited above, many-electron correlations 
are quite essential for the correct description of various characteristics of the cluster 
systems. 

In spite of the fact that the jellium model with all its modifications is rather 
successful in explaining numerous phenomena involving metal clusters it obviously has 
its limits, because this model does not take into account the detailed ionic structure of 
clusters. The correspondence between predictions of the jellium model and the results of 
more advanced quantum chemistry calculations have not been performed in a systematic 
way so far. Partially, this is connected with the fact that quantum chemistry calculations 
are usually limited by small sizes of clusters, while the jellium model becomes adequate 
for larger cluster systems. Knowledge of the ranges of applicability of the jellium model 
and the level of its accuracy is important, because the jellium model often gives much 
more efficient theoretical basis particularly, when dealing with larger cluster systems. 

In this paper we have undertaken detailed systematic theoretical study of structure 
and properties of sodium clusters beyond the jellium model using all-electron ab 
initio theoretical methods based on the Hartree-Fock approximation, density functional 
theory and perturbation theory, for clusters that size is large enough for jellium 
calculations. Namely, we have calculated optimized geometries of neutral and singly- 
charged sodium clusters consisting of up to 20 atoms, their multipole moments (dipole 
and quadrupole), static polarizabilities, binding energies per atom, ionization potentials 
and frequencies of normal vibration modes. We compare results of our calculations with 
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the available experimental data, results of other theoretical works performed both within 
the framework of the jellium model and beyond, using quantum chemistry methods, and 
elucidate the level of accuracy of different theoretical approaches. Also, we demonstrate 
the great role of many-electron correlations in the formation of structure and properties 
of small metal clusters. Our results form a good basis for the detailed study of dynamic 
properties of small metal clusters as well as structure and properties of other atomic 
cluster systems. 

Our calculations elucidate the level of accuracy of various theoretical schemes for 
the treatment of electronic structure in metal clusters, which is important to know and 
is not obvious in advance due to complexity of theoretical methods involved. Some 
characteristics (dipole and quadrupole moments or spectra of normal vibration modes, 
for example), which we have calculated in this paper are new and were not studied 
before, at least according to our knowledge. These characteristics, however, might be, 
important, for instance, when considering dynamics of a cluster beam in an external 
non-homogeneous electric or magnetic field. Indeed, namely, cluster multipole moments 
should be responsible for the cluster isomers separation in the non-homogeneous external 
fields. We analyse the connection between the principal values of the cluster quadrupole 
moments tensor and the cluster shape (oblate, prolate or triaxially deformed). 

The frequencies of the surface and volume vibration modes have been determined 
in the spectra of the cluster normal vibration frequencies and their correspondence to 
the predictions of the dynamical jellium model f4"I|,|4"2f was established. 

Our calculations have been performed with the use of the Gaussian 98 software 
package [5B|. We have used the atomic system of units in this paper, h = m e = |e| = 1 
unless other units are not indicated. 



2. Theoretical methods 



In this work we are studying structure and properties of small sodium clusters on the 
basis of all-electron ab initio many-body theory methods. We calculate the optimized 
geometries of clusters consisting of up to iV < 20 atoms, where N is the number of atoms 
in the cluster. For the sequence of clusters with N < 20, we determine size dependence 
of the cluster ionization potentials, total energies, multipole moments (dipole and 
quadrupole), bonding distances and dipole polarizabilities. We also calculated and 
analyze vibration spectra of the clusters. 

We have done these calculations using different theoretical schemes. We have 
calculated cluster characteristics in the all-electron Hartree-Fock approximation. This 
approximation does not take into account many-electron correlations in the system, 
which turn out to play essential role in the formation of clusters properties. Therefore, 
we also calculate all the characteristics using post Hartree-Fock theories accounting 
for many-electron correlations. Namely, this was done in the M0ller and Plesset 
perturbation theory of the second and the fourth order and the three parameter 
Becke's gradient-corrected exchange functional with the gradient-corrected correlation 
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functional of Lee, Yang and Parr. 

Below, we discuss theoretical methods used in our work. The aim of this discussion 
is to present essential ideas of the methods and give the necessary references, rather 
than to describe them in detail. 



2.1. Hartree-Fock method (HF) 

In the Hartree-Fock approximation, the many-electron wave function of a cluster is 
expressed as antisymmetrized product of the single-electron wave functions, ipi, of cluster 
electrons, which are also often called molecular orbitals. The Hartree-Fock equation for 
the determination of the molecular orbitals ipi reads as (see e.g. 



(-A/2 + U ions + U HF ) A = erfi- (1) 

Here, the first term represents the kinetic energy of the i-th electron, and Ui ons describes 
its attraction to the ions in the cluster. The Hartree-Fock potential Uhf represents the 
Coulomb and the exchange interaction of the electron i with other electrons in the 
cluster, £j is the single electron energy. 

In Gaussian 98, the molecular orbitals, ipi, are approximated by a linear 
combination of a pre-defined set of single-electron functions, \^ known as basis 
functions. This expansion reads as follows: 

N 

^i = ^2 c fiXft, ( 2 ) 

where coefficients c^j are the molecular orbital expansion coefficients, N is the number 
of basis functions, which are chosen to be normalized. 

The basis functions Xn are defined as linear combinations of primitive gaussians: 

v 

where d w are fixed constants within a given basis set, the primitive gaussians, g p = 
g(a,r), are the gaussian-type atomic functions having the following form: 

g(a,r) = cx n y m z l e- ar2 (4) 

Here, c is the normalization constant. The choice of the integers n, m and / defines the 
type of the primitive gaussian function: s, p, d or f (for details see [plfl). 

Substituting these expansions in the Hartree-Fock equations (|]), one can rewrite 
them in the form, known also as the Roothaan and Hall equations: 

N 

Y^i H ^ - £iSp,)cvi = n=l,2,...,N (5) 

Being written in the matrix form, this equation reads as: 

HC = SCe, (6) 
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where each element is a matrix. Here, e is a diagonal matrix of orbital energies, each of its 
elements £j is the single-electron energy of the molecular orbital ipi, H is the Hamiltonian 
matrix as it follows from (01), S is the overlap matrix, describing the overlap between 
orbitals. For more details regarding this formalism see . 

Equations (0) are none linear and must be solved iteratively. The procedure which 
does so is called the Self- Consistent Field (SCF) method. 

The above written equations consider the restricted Hartree-Fock method. For the 
open shell systems, the unrestricted Hartree-Fock method has to be used. In this case, 
the alpha and beta electrons with spins up and down are assigned to different orbitals, 
resulting in two sets of molecular orbital expansion coefficients: 

N 
N 

# = (7) 

The two sets of coefficients result in two sets of the Hamiltonian matrices and the 
two sets of orbitals. 

2.2. M0ller-Plesset perturbation theory method (MP n ) 

The Hartree-Fock theory provides an inadequate treatment of electrons motion within 
a molecular system, because it does not properly treat many-electron correlations. The 
many electron correlations can be accounted for using different methods. The most 
straightforward way for achieving this goal is based on the perturbation theory. Indeed, 
the total Hamiltonian, H, of the cluster can be divided into two parts 

H = H + V (8) 

Here H is the Hamiltonian corresponding to the Hartree-Fock level of theory and V is 
the residual interelectron interaction, which can be treated as a small perturbation. 

Considering V as a small perturbation one can construct the solution of the 
Schrodinger equation for many-electron system in an arbitrary order of the perturbation 
theory. The perturbation theory of this type is well known since the work by M0ller- 
Plesset []59| and can be found in numerous textbooks on quantum mechanics (see 



e 



•g- m 



Below we refer to this theoretical method as to the M0ller-Plesset perturbation 
theory |59| of the second or forth order, MP 2 or MP 4 . Indices here indicate the order 



of the perturbation theory. 
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2.3. Density functional methods (B3LYP) 

The density functional theory (DFT) is based upon a strategy of modelling electron 
correlation via general functionals of the electron density. Within the DFT one has to 
solve the Kohn-Sham equations, which read as (see e.g. p|,P,p!THl4|) 

where the first term represents the kinetic energy of the i-th electron, and Ui ons describes 
its attraction to the ions in the cluster, Vh is the Hartree part of the interelectronic 
interaction: 

V H {r)= [ j^ljrdr', (10) 
and p{f ') is the electron density: 

N 

where V xc is the local exchange-correlation potential, ipi are the electronic orbitals and 
N is the number of electrons in the cluster. 

The exchange-correlation potential is defined as the functional derivative of the 
exchange- correlation energy functional: 

V. = ^M, (12) 

The approximate functionals employed by DFT methods partition the exchange- 
correlation energy into two parts, referred to as exchange and correlation parts: 

E xc [p] = E x (p) + E c (p) (13) 

Physically, these two terms correspond to same-spin and mixed-spin interactions, 
respectively. Both parts are the functionals of the electron density, which can be of 
two distinct types: either local functional depending on only the electron density p or 
gradient- corrected functionals depending on both p and its gradient, Vp. 

In literature, there is a variety of exchange correlation functionals. Below, we refer 
only to those, which are related to the calculation performed in this work. 

The local exchange functional is virtually always defined as follows: 

E' DA = -\^ f p^dh (14) 

This form was developed to reproduce the exchange energy of a uniform electron 
gas. By itself, however, it is not sufficient for the adequate description of atomic clusters. 

The gradient-corrected exchange functional introduced by Becke |JT] and based on 
the LDA exchange functional reads as: 

^LDA _ / P X 



E m* = e lda _ / ; " d\ (15) 
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where x=p" 4 / 3 |Vp| and 7 = 0.0042 is a parameter chosen to fit the known exchange 
energies of the noble gas atoms. 

Analogously to the above written exchange functionals, there are local and gradient- 
corrected correlation functionals, for example, those introduced by Perdew and Wang 
|62j or by Lee, Yang and Parr [|63|. Their explicit expressions are somewhat lengthy and 



thus we do not present them here and refer to the original papers. 

In the pure DFT, an exchange functional usually pairs with a correlation functional. 
For example, the well-known BLYP functional pairs Becke's gradient-corrected exchange 
functional with the gradient-corrected correlation functional of Lee, Yang and 
Parr |53|. 



In spite of the success of the pure DFT theory in many cases, one has to admit 
that the Hartree-Fock theory accounts for the electron exchange the most naturally 



and precisely. Thus, Becke has suggested [61] functionals which include a mixture of 



Hartree-Fock and DFT exchange along with DFT correlations, conceptually defining 

E xc as! 



= chfB%* + c DFT E£ \ (16) 

where chf an d cdft are constants. Following this idea, a Becke-type three parameter 
functional can be defined as follows: 



E 



B3LYP 



E L X DA + c (E» F - E L X DA ) + c x (E* 88 - E L X DA ) + 



E 



VWN3 



+ cJE, 



LYP 



£VWN3\ 



(17) 



Here, Co = 0.2, c x = 0.72 and c c = 0.81 are constants, which were defined by 
fitting to the atomization energies, ionization potentials, proton affinities and first-row 



atomic energies f58j. E^ DA and E^ 88 are defined in ([14]) and fll5|) respectively. E^ F is 



the functional corresponding to Hartree-Fock equations ([!]). The explicit form for the 
correlation functional E VWN3 as well as for gradient-corrected correlation functional of 
Lee, Yang and Parr, E^ YP , one can find in [64]] and [63| correspondingly. Note that 



instead of E^f WN3 and E^ YP in ( |P7D one can also use the Perdew and Wang correlation 



functional [62]. 



2.4- Geometry optimization 

The cluster geometries, which we have calculated in our work, have been determined 
using the geometry optimization procedure. This procedure implies the calculation 
of the multidimensional potential energy surface for a cluster and then finding local 
minima on this surface. The key point for this search is fixing the starting geometry of 
the cluster, which could converge during the calculation to the local or global minimum. 
There is no unique way in achieving this goal with Gaussian 98. 

In our calculations, we have created the starting geometries empirically, often 
assuming certain cluster symmetries. Note, that during the optimization process the 
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geometry of the cluster as well as its initial symmetry sometimes change dramatically. 
All the characteristics of clusters, which we have calculated and presented in next 
section, are obtained for the clusters with optimized geometry. 

In our calculations, we have made no assumptions on the core electrons in the 
optimized clusters, which means that all electrons available in the system, have been 
taken into account, when computing potential energy surface. For clusters with N > 10, 
this process becomes rather computer time demanding. Thus, in this work we have 
limited our calculations by clusters consisting up to N < 20. 



2.5. Normal vibrations 

Knowledge of the potential energy surface in the vicinity of a local minimum, allows 
one easily to determine corresponding normal vibration modes of the system. We have 
performed such calculation and determined the vibration energy spectrum for a number 
of clusters. Particular attention in this calculation has been paid to the identification 
of the breathing and the surface vibration modes and comparison their frequencies with 
those predicted in [IT, 12" for spherical sodium clusters on the basis of the dynamical 
jellium model. 



3. Results of calculations and discussion 



In this section we present the results of calculations performed with the use of 
methods described above. We have calculated the optimized geometries of neutral and 
singly charged sodium clusters consisting of up to 20 atoms, their multipole moments 
(dipole and quadrupole), static polarizabilities, binding energies per atom, ionization 
potentials and frequencies of the normal vibration modes. We compare results of our 
calculations with the available experimental data and the results of other theoretical 
works performed both within the framework of the jellium model and beyond, using 
quantum chemistry methods and establish the level of accuracy of different theoretical 
approaches. Particular attention is paid to the clusters in the range 10 < N < 20, 
because some characteristics of the clusters in this size range have been calculated on 
the ab initio basis in our paper for the first time. Also, we demonstrate the great role of 
many-electron correlations in the formation of structure and properties of small metal 
clusters. 



3.1. Geometry optimization of Na n and Na£ clusters 

Results of the cluster geometry optimization for neutral and singly charged sodium 
clusters consisting of up to 20 atoms shown in figures |TJ and ^] respectively. The cluster 
geometries have been determined using the methodology described in section |2|. Namely, 
the optimization of the cluster geometries has been performed with the use of B3LYP 
and MP 2 methods. 
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Parte. 



Figure 1. Optimized geometries of neutral sodium clusters Na 2 — Naio (part a), 
Nan — Nais (part b) and Naig — Na 2 o (part c). The interatomic distances are given 
in angstroms. The label above each cluster image indicates its point symmetry group 
and the calculation method by which the cluster was optimized. 
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Figure 2. Optimized geometries of singly charged sodium clusters Na 2 — Na^ (part 
a) and Na^ 2 — Na 21 (part b). The interatomic distances are given in angstroms. The 
label above each cluster image indicates the point symmetry group and the calculation 
method by which the cluster was optimized. 
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For clusters with N < 6, we preferably used the MP 2 method. This method leads 
to the results, which are in a reasonable agreement with those derived by other methods 
(see e.g. [p!j , pjl ). For example, the side bond length in the rhomboidal Na^ cluster 
calculated in |]T6[ by the all-electron Hartree-Fock method is equal to 3.74 A, while in 
our case it is equal to 3.56 A. The smaller diagonal value for Na^ is equal to 3.25 A 
in 



TB|, while we determine it as 3.18 A. 



The MP2 method becomes more and more computer time demanding with the 
growth cluster size. This happens due to increase in a number of integrals involved in 
the computations. It turns out that for larger cluster systems the B3LYP method is 
more efficient. The accuracy of the B3LYP method is comparable with the accuracy 
of the MP2 method, as it is clear from the comparison of the B3LYP and MP 2 cluster 
geometries with those computed in |16| by the configuration interaction method. 

Clusters of a certain size can possess various isomer forms, those number grows 
dramatically with increasing cluster size. We illustrate the situation, and calculate 
several isomers of the the Na^, Na^, Na w , Na n and Na 2 o clusters. They all are 
presented in figure [I]. Note, that the linear and equilateral triangular Na% isomers, 



have not been described in the earlier papers fl5|-4T7| (see also [|Tl], |13|,|14[), in which 
isosceles triangular isomers were considered. The comparison of properties (dipole and 
quadrupole moments, total energies, bonding distances) of these clusters will be given 
below. 

On the example of the Na^ cluster, we demonstrate how the multiplicity of an 
electronic state of the system can influence its geometry. Figure |I] shows that the Na^ 
cluster has the rhomboidal geometry corresponding to the D 2 h point symmetry group, 
if the multiplicity of the cluster is equal to 1, while, for the multiplicity being equal to 3, 
the cluster has the quadratic geometry characterised by the point symmetry group. 

Sodium clusters with N < 5 have the plane structure, while for N = Q both plane 
and spatial isomers are possible. This feature is consistent with the jellium picture and 
can be explained from the minimization principle for the cluster surface. Indeed, the 
surface of small plane cluster isomers is less in comparison with the surface of their 
possible spatial forms. 

Comparison of geometries of the neutral and singly- charged clusters presented in 
figures [l] and |2| shows their significant difference. For smaller sizes (N < 8), singly- 
charged and neutral clusters have sometimes different point symmetry groups and 
bonding distances (see images of the Na^, Na$, Na^ and Na$ clusters and their ions). 
The alteration in the geometry of cluster ions occurs due to the excessive positive charge 
available in the system. The structural change of cluster ions becomes less profound with 
increasing cluster size, see clusters with N > 10, because the excessive positive charge in 
this case turns out to be insufficient to produce substantial change in a massive cluster, 
although sometimes (compare iVais and iVa^) noticeable change in the cluster geometry 
is also possible. 

The striking difference in geometries of small singly charged and neutral clusters 
is closely linked to the problem of cluster fission. It is natural to assume that with 
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increasing cluster charge small clusters should become unstable and fragment into two 
parts, while for larger cluster sizes one can expect quasi-stable configurations, which 
should decay via the fission process. Calculation of such configurations is an interesting 
task, because it may provide the essential information on the predominant fission 
channels in the system. We do not perform such an analysis in our work, but draw 
attention that geometries of the cluster ions, like Na^ , Na£, Na$ and Na± 5 , lead to 
the obvious hints on the possible fragmentation channels in these cluster systems. 

Figure [l] shows that the clusters Nag and Na2o have the higher point symmetry 
group Tfi as compared to the other clusters. This result is in a qualitative agreement with 
the jellium model. According to the jellium model |36|-^0|, clusters with closed shells 
of delocalized electrons have the spherical shape, while clusters with opened electron 
shells are deformed. The jellium model predicts spherical shapes for the clusters with 
the magic numbers N = 8, 20, 34, 40..., having respectively the following electronic shells 
filled: ls 2 lp 6 , ld 10 2s 2 , l/ 14 , 2p 6 , 

We have also found the Td symmetry group isomer for the Na\o cluster. However, 
this cluster isomer is not the lowest energy isomer of Na w (see table [SI]). The similar 
situation occurs in the jellium model, where the ls 2 lp 6 , 2s 2 closed shell electronic 
configuration does not minimize the cluster total energy. 

Note also, that both the LDA and HF jellium models predict some deviation from 
sphericity for the Na\s cluster [|(J having Id subshell filled, which is a result of electron 
configurations mixing. This fact is also in a qualitative agreement with the results of 
our ab initio calculations. The point group symmetry of the Nai$ cluster, C$ v , is lower 
than T d , which is the point symmetry group for the Na 8 and Na 2 o clusters, and even 
lower than the point symmetry group for some opened shell clusters, like Na? and Nai$, 
having the point symmetry group D^. 

Note that there are some clusters possessing relatively low point symmetry group, 
that nevertheless is quite close to the higher point symmetry group. The higher 



symmetry breaking is not occasional and can be explained via the Jahn- Teller effect |60 
Such situation occurs, for example, in the Na 9 and Nan clusters, which posses the C 2v 
point symmetry group, but their geometry is close to the geometry of the D 3h group. 

The jellium prediction on the sphericity of the magic clusters works not so well 
for cluster ions. Indeed, the geometry and the point symmetry group of Na$ does not 
allow one to state the higher sphericity of this cluster as compared to its neighbours. 
The analysis of the quadrupole moments and cluster deformations performed below 
demonstrates this conclusion quite clearly. This happens because forces emerging in 
the cluster during its transition from neutral to singly charged state turns out to be 
insufficient to rearrange the cluster geometry from deformed to spherical one. 

We have found two isomers of the Na2o cluster, which have rather regular structure 
and differ significantly one from another. The cluster geometries presented in figure 
[TJ allow one to assume that there exist at least two independent paths of the cluster 
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Figure 3. Averaged bonding distance as a function of cluster size for optimized 
geometries of neutral sodium clusters. For some cluster numbers more than one isomer 
has been considered. In these cases, labels indicate the point symmetry group of the 
corresponding isomers. Geometries of the optimized clusters one can find in figure |l]. 

structure formation. Indeed, the following isomers 

Na$ 5v -> Na 7 -> Na^ v -> Na 13 -> Na 15 -> Na w -> Na 17 -> Na 18 -> Na 19 -> Na^ v 

probably belong to the chain leading to the formation of the C 2v isomer of the Na2o 
cluster, while the clusters 

Na° 3h -> Na 8 -> Na 9 -> JVog -> iVag -> iVa 12 -> iVai 4 -> iVa^ 

form the path on which the isomer of the Na 2 o cluster is formed. Figure |I] clearly 
shows the steps of the cluster formation process along these two paths. Although, for 
most of N, we have calculated isomers belonging to one path or another, it is natural 
to assume that the two different type of geometries exist for all N, similar to how 
it happens for Na$ and Na2o clusters. For clusters smaller than Na^, one can not 
distinguish the two paths clearly enough as it is seen from figure 0. Conclusions made 
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Figure 4. Averaged bonding distance as a function of cluster size for optimized 
geometries of singly charged sodium clusters. For some cluster numbers more than 
one isomer has been considered. In these cases, labels indicate the point symmetry 
group of the corresponding isomers. Geometries of the optimized clusters one can find 
in figure ^. 

for neutral clusters regarding the growing process are applicable to the great extent to 
singly charged cluster ions as it is clear from figure |2|, although cluster ions geometries 
sometimes differ substantially from their neutral prototypes. 

Cluster geometries allow one easily to compute and analyze the average bonding 
distance as a function of cluster size. The result of this analysis for neutral and singly 
charged sodium clusters is presented in figures |3| and [|. These figures show how the 
average bonding distance converge to the bulk limit indicated in the figures by horizontal 
lines. When calculating the average bonding distance in a cluster, interatomic distances 
smaller than 4.1 A have only been considered. This upper limit on the interatomic 
distances has been chosen as a distance, which is 10 per cent larger than to the bcc- 
lattice nearest neighbour distance in the bulk sodium. 

Figures |3| and [| show that the dependence of the average bonding distance, (R) , on 
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Figure 5. The principal values of tensor i£y for optimized neutral sodium clusters 
as a function of cluster size calculated by the B3LYP method. Squares, circles and 
triangles represent the R xx , Ryy and R zz tensor principal values respectively. For some 
clusters, more than one isomer has been considered. In these cases, labels indicate the 
point symmetry group of the corresponding isomers. Geometries of the optimized 
clusters one can find in figure || 



cluster size is non-monotonous. For neutral clusters, one can see odd-even oscillations 
of (R) atop its systematic growth and approaching the bulk limit. These features 
have the quantum origin and can be explained by the derealization of valence atomic 
electrons. Indeed, the odd-even oscillations arise due to the spin paring of the delocalised 
electrons. This type of behaviour is also typical for other cluster characteristics and will 
be discussed below in more detail. Relatively large increase of the average distance, seen 
for small sodium cluster ions with N < 9, is also qualitatively clear. It can be explained 
by the Coulomb instability developing in the cluster with increasing its ionization rate. 

Cluster shape can be characterized by the oblate, prolate or triaxial deformation. 
The prolate deformation of the cluster is characterized by larger distortion of the ionic 
charge distribution along z-axis as compared to distortions along x- and y axes. In 
the oblate deformation case the situation is opposite. Deformations of the ionic charge 
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Figure 6. The principal values of tensor i?y for optimized singly charged sodium 
clusters as a function of cluster size calculated by the B3LYP method. Squares, circles 
and triangles represent the R X x, Ryy and R zz tensor principal values respectively. For 
some clusters, more than one isomer has been considered. In these cases, labels indicate 
the point symmetry group of the corresponding isomers. Geometries of the optimized 
clusters one can find in figure |[ 



distribution in x- and y- directions are larger than in z-direction. In both cases the 
deformations along x- and y- directions are equal. The triaxial shape deformation is 
characterized by unequal distortions of the ionic charge distribution along x-, y- and z- 
directions. Often, however, two of three deformations are close to each other and this 
allows one to discuss the triaxially deformed prolate or oblate cases. Knowledge of the 
type of the cluster deformation is quite useful for the comparison with the jellium model 
results and the analysis of the metal cluster photon absorption spectra by metal clusters 
(see 0). 

The type of cluster deformation can be easily determined by the principle values of 
the tensor R^j = Y2 x i x j- Here, the summation is performed over all ions in the system. 
The principle values of this tensor R xx , R yy and R zz define the dimensions R x , R y and 
R z of the ionic charge distribution in the cluster along the principle axes x, y and z via 



Structure and properties of small sodium clusters. 



21 



the relations: R x = ^JR XX /N \ R y = ^R yy /N and R z = a/ R zz /N. Note that tensor 
Rij is closely connected with the cluster moment of inertia tensor and the quadrupole 
moment tensor of the ionic distribution. 

In figures [|and [5] we present the principle values R xx , R yy and R zz for a sequence of 
neutral and singly charged clusters respectively Figures | and [] demonstrate how the 
cluster deformation change as a function of cluster size. Figure |5] shows that all three 
principle values are equal for the tetrahedron group isomers of the magic clusters Na 8 
and Na2o- This feature is in the qualitative agreement with the jellium model, which 
predicts spherical shapes for the magic clusters. In many cases two of three principal 
values of R^ are equal or nearly equal. Using the definition of the prolate and oblate 
cluster distortions given above and figures |5] and || one can easily determine the type of 
cluster deformation. For example, clusters Na2, Naio, Nai$ and iVaig have the prolate 
deformation along z-principle axis, because the following condition R xx = R yy < R zz is 
fulfilled. The clusters Na 6 and Na 7 possess the prolate deformation because in this case 
Rxx = Ryy > Rzz- Figures H and | show that most of clusters are triaxially deformed. 
However, it is often possible to assign clusters triaxially deformed prolate or oblate 
shape, because two of three principle values are close to each other. Thus, for instance, 
Na4, Nciis are triaxialy prolate clusters, while Na u is a triaxialy oblate one. Figures [| 
and |6] also show the relative value of prolate and oblate deformations in various clusters. 

One can define a tensor analogous to R^, but for electrons. We do not plot the 
principle values of such a tensor because they are very close in absolute value to the 
principle values shown in figures [5| and |6] and could be traced from the principle values 
of the cluster total quadrupole moment tensor considered below in subsection |3.4| . 

3.2. Binding energy per atom for small neutral and singly-charged sodium clusters. 

The binding energy per atom for small neutral and singly- charged sodium clusters is 
defined as follows: 



E b /N = E l - E N /N (18) 
E+/N = ((N - 1)^ + Ef - E+) /N, (19) 

where E^ and E^ are the energies of a neutral and singly- charged N-atomic cluster 
respectively. E\ and E± are the energies of a single sodium atom and an ion. 

Figures [7| and || show the dependence of the binding energy per atom for neutral and 
singly-charged clusters as a function of cluster size. The energies of clusters have been 
computed using the B3LYP, MP^ and HE methods described in section |2|. For clusters 
with N < 8, computations of the energies have been performed by the three methods 
for the sake of comparison. We wanted to compare the methods by their accuracy and 
computation efficiency. The results of our calculations have also been compared with 
those derived by the configuration interaction (CI) method in [|l5| -|T7|j). Figures |7| and 
|8| demonstrate that the results of the MP4 and B3LYP methods are in a reasonable 
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Figure 7. Binding energy per atom for neutral sodium clusters as a function of 
cluster size. Circles represent the binding energies per atom calculated by the B3LYP 
method, lower and upper triangles correspond to the energies obtained by the MP4 
method and in the HF approximation respectively. Squares show the result of the 
configuration interaction approach from the work by Bonacic-Kotecky et al (for details 
see fL6| , |lg[| ), Some points in figure have labels, indicating the point symmetry group 
of the isomers represented. Geometries of the corresponding clusters one can find in 
figure [1} 



agreement with each other and with the CI results. The HF points significantly differ 
from the MP 4 , B3LYP and CI ones, which demonstrates the importance of many- 
electron correlations, taken into account in the MP 4 , B3LYP and CI methods and 
omitted in the HF approximation. Note that the energy of Na2, if computed in the 
pure HF approximation, is close to zero, which means that bonding in this molecule 
takes place mainly due to many-electron correlations. 

The energies of clusters larger than Na$ and Na$ have been computed by the 
B3LYP method only, because this method is more efficient than MP 4 and the accuracy 
of both methods is comparable. 

Figures [7| and |8] demonstrate the even-odd oscillation behaviour in the dependence 
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Figure 8. Binding energy per atom for singly charged sodium clusters as a function of 
cluster size. Circles represent the binding energies per atom calculated by the B3LYP 
method, lower and upper triangles correspond to the energies obtained by the MP4 
method and in the HF approximation respectively. Squares show the result of the 
configuration interaction approach from the work by Bonacic-Kotecky et al (for details 
see Jl6|,[l|]). Some points in figure have labels, indicating the point symmetry group 
of the isomers represented. Geometries of the corresponding clusters one can find in 
figure [2| 

of binding energy on cluster size. Indeed, for singly charged clusters, odd numbers 
corresponding to the singlet multiplicity have higher energies as compared to their even 
neighbours. Analogous situation takes place for neutral clusters. In this case, even 
cluster numbers have higher binding energies as compared to their odd neighbours. Note 
that for neutral clusters this phenomenon occurs simultaneously with slight systematic 
growth of the binding energies per atom with increasing cluster size. 

Figures [7] and [| also show that the binging energy per atom in the magic neutral 
clusters, Na 8 and Na 2 o, is a little higher as compared to other clusters of the close size. 
The similar situation takes place for the Nag cluster in the ionic case. This feature 
can be qualitatively understood on the basis of the jellium model: increasing the magic 
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clusters binding energy takes place due to the delocalised electrons shell closure. Note 
that the binding energy per atom for the magic Na^ turns out to be smaller than 
that for the neighbouring cluster ions. This happens because this particular cluster ion 
isomer is characterized by the Oh point symmetry group. Cluster isomers based on this 
point symmetry group usually have the lower binding energy per atom as compared to 
the isomers based on the icosahedron point symmetry group like those with iV > 13 
shown in figures [l] and |[ 

Tables |A1| and |A2| given in |Appendix A| provide the accurate values of the cluster 



total energies calculated by MP 4 , B3LYP and HF methods. For neutral clusters with 
iV < 8, we also present the cluster energies calculated in |l6j by the CI method. The 
values given in these tables have been used to plot figures [7] and |[ For some clusters, 
energies of different symmetry isomers are also given in the tables. 

3.3. Ionization potentials 

Let us now consider how the ionization potentials of sodium clusters evolve with 
increasing cluster size. Experimentally, such a dependence has been measured for sodium 
clusters in f|^0|. 



The ionization potential of a cluster consisting of N atoms is defined as a difference 
between the energy of the singly- charged and neutral clusters: 

IP = E+- E N (20) 

Figure § shows the dependence of the clusters ionization potential on N. Figure || 
demonstrates the comparison of the results derived by different methods, B3LYP, MP4 



and HF (see section 0), with the experimental data from |8j and |p0|| . The results of 
the B3LYP and MP4 methods are in a reasonable agreement with the experimental 
data, while the ionization potentials calculated on the basis of the HF approximation 
differ substantially from the experimental observations. This comparison shows the role 
of many-electron correlations in the formation of the cluster ionization potentials. The 
correlation effects are taken into account by the B3LYP and MP 4 methods and omitted 
in the HF approximation. 

Figure |9] demonstrates that the ionisation potentials drop with increasing cluster 
size, which is consistent with predictions of the classical spherical droplet model. 
However, this process has many irregularities, which have quantum origin. Indeed, 
the dependencies derived by the MP4 and B3LYP methods as well as the experimental 
one have a prominent odd-even oscillatory tendency. The maxima in these dependences 
correspond to the even-N-clusters, which means their higher stability as compared to the 
neighbouring odd-N-clusters. This happens because the multiplicities of the even- and 
odd-N-clusters are different, being equal to one and two correspondingly. Interestingly 
enough that the B3LYP method reproduces correctly even the experimentally observed 
irregularity in the odd-even oscillatory behaviour, which happens at iV = 16 and N — 17, 
and some other minor details of the experimental data. 
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Figure 9. Ionization potentials of neutral sodium clusters as a function of cluster 
size. Circles show the results derived by the B3LYP method. Triangles and 
rhomboids represent the ionization potentials calculated by the HF and MP4 methods 
respectively. Filled and open squares are the experimental values taken from |^] and || 
respectively. For some clusters, more than one neutral and/or singly charged cluster 
isomer has been considered. In these cases, labels indicate the point symmetry group 
of the initial neutal and the final charged cluster isomers used for the calculation of 
the ionization potential. 



A significant step-like decrease in the ionization potential value happens at the 
transition from the dimer to the trimer cluster and also in the transition from Na$ to 
Na 9 . Such an irregular behaviour can be explained by the closure of the electronic ls- 
and lp-shells of the delocalized electrons in the clusters iVa 2 and Na$ respectively. The 
next significant drop in the ionization potential value takes place in the transition from 
the magic Na2o to the Na2i cluster. 
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Figure 10. Dipole moments of the optimized neutral sodium clusters as a function 
of cluster size calculated by the B3LYP method. For some clusters, more than one 
isomer has been considered. In these cases, labels indicate the point symmetry group 
of corresponding isomers. Geometries of the optimized clusters one can find in figure 
0. 1 Debye=0.3935 a.u. 



S.Jf.. Multipole moments 

We have calculated multipole moments (dipole, quadrupole, octapole and hexadecapole) 
for the sodium clusters those geometry is shown in figures [l] and |^. In figures [H] and [Ll], 
we plot the absolute values of the dipole moments for the neutral and singly charged 
sodium clusters clS db function of cluster size. 

The dipole moments of some sodium clusters (see figure [10]), which we predict in our 
paper, arise due to the fact that the electron charge distribution not always matches the 
ionic charge distribution and can be shifted with respect to the cluster centre of mass. 
Our calculations show that only clusters with the C-point symmetry groups, like the 
isosceles triangle isomers of Na 3 , the pentagonal Na 6 pyramid isomer, Nai 2 , iVa 18 and 
others, possess dipole moments. These clusters have either an axis of a certain order or 
a plane of symmetry, but no perpendicular symmetry elements (plains or axes). This 
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Figure 11. Dipole moments for the optimized singly charged sodium clusters as a 
function of cluster size calculated by the B3LYP method. For some clusters, more 
than one isomer has been considered. In these cases, labels indicate the point symmetry 
group of corresponding isomers. Geometries of the optimized clusters one can find in 
figure |2| 



rule remains correct even for the Na2o cluster isomer with the symmetry Ci v , which 
has the closed shell configuration ls 2 lp 6 ld 10 2s 2 of delocalised electrons according to the 
jellium model. Geometries of the cluster ions differ significantly from the geometries of 
the corresponding neutral clusters, but the rule formulated above on the appearance of 
the cluster dipole moments remain valid in this case also as it is clear from figure |ll|. 
The principal values of the quadrupole moments tensor for the optimized neutral 



and singly charged clusters are presented in figures [12] and [13| respectively. For clusters 
with an axis of symmetry, this axis has been chosen as z-axis of the coordinate system, in 
which the calculation of the quadrupole moments has been performed. The quadrupole 
moment tensor is defined as an average value of the following operator: 

Qij = q(3xiXj - 5ijr 2 ) (21) 

Here, the summation is performed over all electronic and ionic charges in the cluster. 
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Figure 12. The principal values of quadrupole moment tensor for the optimized 
neutral sodium clusters as a function of cluster size calculated by the B3LYP method. 
Squares, circles and triangles represent the Q xx , Qyy and Q zz tensor principal values 
respectively. For some clusters, more than one isomer has been considered. In these 
cases, labels indicate the point symmetry group of corresponding isomers. Geometries 
of the optimized clusters one can find in figure |l|. 

Note that the trace of the tensor is equal to zero. 

The ionic part of Qij can be expressed via the components of the tensor Rij discussed 
in section Note that the knowledge of and allows one to construct easily 
the tensor analogous to R^, but for electrons. This might be useful for the analysis of 
deformations of electron density distribution in a cluster. 

The quadrupole moment tensor can be expressed via the tensor = (J2Q. x i x j)i 
characterising the averaged dimensions of the total charge distribution. Here, brackets 
mean averaging over the electronic charge distribution. The principal values of the tensor 
Qij should be negative at least for neutral clusters, because electron density is spilled 
out of the cluster, which makes its distribution a little broader than the distribution of 
ions. The similar situation takes place for cluster ions, but in this case there is non- 
compensated positive charge in the system, which brings certain positive contribution 
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Figure 13. The principal values of quadrupole moment tensor for the optimized singly 
charged sodium clusters as a function of cluster size calculated by the B3LYP method. 
Squares, circles and triangles represent the Q XXl Q yy and Q zz tensor principal values 
respectively. For some clusters, more than one isomer has been considered. In these 
cases, labels indicate the point symmetry group of corresponding isomers. Geometries 
of the optimized clusters one can find in figure |^. 

to Qij and makes the principal values of Qij positive in some cases. 

The numerical analysis performed in this work shows that for neutral sodium 
clusters the principal values of Qij are always negative, while for the small cluster ions: 
Na>2, Naf and Na^ (C^), some of the principal values are positive. 

The principle values of the quadrupole moment tensor characterize the distortion 
of the total cluster charge distribution. Indeed, figure [L^ shows that the Na$ and 
Na 2 Q tetrahedron group isomers have the zero quadrupole moments, which reflect the 
closeness to sphericity of the magic clusters. Our calculations demonstrate that for 
some open shell clusters like Nan and Nau the quadrupole moments turn out to be 
rather small, although the ionic charge distribution in these clusters has the prominent 
deformation as it is clear from figures [I] and |5|. The small quadrupole moments in these 
clusters is the result of compensation of the electron and ion components of Qij . 
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The quadrupole moments diagram allows one to make some conclusions on the 
type of the shape of the total charge distribution in a cluster. The averaged dimensions 
of the cluster total charge distribution in x-, y- and z- directions can be characterized 
by quantities Q\ = Q zz = (J2 ez 2 }, Q^ = Q xx = (J^ex 2 ) and = Q yy = (X>2/ 2 )- 
Here, the summation is performed over all electrons and ions in the cluster and brackets 
mean averaging. These quantities are connected with the quadrupole moments tensor 
defined in (pi]). Indeed, in both the prolate and oblate cases, when = Qy = Q 1 - and 
Qz = QK the principal values of the tensor Qij read as 



Qzz = 2(Q" - Q x ) 

Q XX = (Q ± -Q ll ) = -^ 

Qyy = Q XX = "% (22) 

These equations define the important relationships between the principal values of 
the quadrupole moments tensor in the oblate and prolate cases and help understanding 
the quadrupole moments diagrams shown in figures ^ and |I3[ 

Equations (|22| ) show that the sign of the principal values Q xx , Q yy and Q zz depends 
on the relative value of Q an d Q L - With the use of equations (^) and the cluster 
quadrupole moment diagrams shown in figures [12] and [13|, one can easily analyse the 
total charge distribution of the clusters shown in figures [I] and |2|. Note that conclusions 
made on the shape of the total charge distribution and the shape of ionic component 
(see figures |5| and |6|) sometimes differ significantly one from another for some clusters. 
For example, the ionic charge distribution in the Na\2 cluster has a prolate shape, while 
the total charge distribution is oblate. 

The quadrupole moments of singly charged sodium clusters differ substantially from 
those for the corresponding neutral ones. The excessive positive charge leads to the 
rearrangement of the cluster structure and to the appearance of the quadrupole moment 
in the cluster ions like Na% and iVa^ . Although, the electron exchange-correlation force 
in a cluster turns out to be insufficient to change the cluster geometry so significantly 
to make the magic cluster ion Nag, having the closed shell electronic structure of 
delocalised electron, spherical-like without quadrupole moment. Instead, Nag remains 
a noticeable deformation. 

Let us now discuss the idea for which the cluster multipole moments play the 
crucial role and consider the possibility of the cluster isomers separation by placing 
the mass selected cluster beam in the inhomogeneous external field. As we have seen 
from the calculations presented above, different cluster isomers of the same mass often 
possess different structure and as a result of that different multipole moments (dipole 
or quadrupole). However, such cluster isomers are indistinguishable in the nowadays 
experiments with mass selected cluster beams. They can nevertheless be separated if 
one puts the mass selected cluster beam in the inhomogeneous external field. Let us 
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estimate this effect for the characteristic values of the dipole and quadrupole moments 
calculated above. 



From the dipole moments diagrams shown in figures 10 and 11 one can conclude 
that the difference in dipole moments for some cluster isomers can be as large as IDebye 
and for the quadrupole often it is about AODebye-X or even larger. The force acting on 
the cluster with the dipole moment D in an external inhomogeneous electric field E(r) 
is equal to f65j 

F D (r) = V{D • E(r)}. (23) 

The components of the force acting on the cluster with quadrupole moment Qij is as 
follows [05(1 



if (r) = V^V^r)}. (24) 

Here, the summation is assumed over the repeated indices j and k of the vector and 
tensor components in the right hand side of (|24]). 

Let us introduce the time period r during which the cluster beam passes the 
inhomogeneous electric field. One can estimate the distance A on which isomers will be 
separated during this period of time as A ~ Ft 2 /2M, where M is the mass of the isomer 
considered and F is the force acting on either the dipole (see (|23|) ) or quadrupole (see 
(H)) moment of the cluster. Substituting in these equations the characteristic values 
for the dipole and quadrupole moments, assuming that the inhomogeneity of the electric 
field is about VE ~ 5 • 10 3 V/cm 2 , one derives from (^) (|24j) that during the period 

10" 3 s the isomers with N = 3 and 5D ~ IDebye become separated on A ~ 0.7mm 
and that A ~ 2.8mm for SQ ~ AODebye-A, r ~ 10s, N = 5 and no dipole moment. 

These estimates demonstrate that one can create significant separation distances 
for reasonably short periods of time with the electric field strengths and their gradients 
achievable in laboratory conditions. The experiments with mass selected and isomer 
separated cluster beams could provide the most accurate information on the structure 
and properties of atomic clusters. 

3.5. Polarizabilities 

We have calculated the polarizabilities for the optimized neutral sodium clusters (see 
figure |l|) as a function of cluster size. Results of this calculation are shown in figure [14|. In 



this figure, we also plot experimental points from [21]. Calculation of the polarizabilities 
has been performed by the B3LYP method. Figure |14] demonstrates quite reasonable 
agreement of the B3LYP results with the experimental data. 

In figure we also compare the polarizabilities calculated in our work with 
those derived by other theoretical methods [19|, |32| . This figure demonstrates a 
satisfactory agreement of the results of different approaches with each other and with 
the experimental data. This comparison is quite important, because in our work as well 
|19| the polarizabilities have been calculated using all electron ab initio approach, 



as m 



while in p2j they were obtained with the use of pseudopotentials. Note that our points 
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Figure 14. Static mean polarizability per atom for neutral sodium clusters normalized 
to the polarizability of a single sodium atom. Circles show the results derived in this 
work by the B3LYP method. For some clusters, more than one isomer has been 
considered. In these cases, labels indicate the point symmetry group of corresponding 
isomers. Stars and triangles represent the polarizabilities calculated in and p2| ] 
respectively. Squares are the experimental values taken from |2lf| . 



are closer to the experimental values than those from [[B|, in spite of the fact that 
both calculations have been performed on the basis of the density functional theory. 
The difference between the two schemes of calculation arise in the form of the density 
functional and the emploied set of the basis functions. In [TT9[| , the so-called Perdew- 
Wang-91 density functional |62|] was used, while we applied its B3LYP form. 

Let us also compare the polarizabilities for the Na$ and Na2o clusters calculated 
in the random phase approximation with exchange in the spherical jellium model, 
c*Na$ = 755a. u. and aNa 20 — 1808a.it. [f)6|| , with our results: «Ar ag = 797 'a.u. and 
c*Na 20 = 1964a. u. The closeness of the values show that the detailed ionic core structure 
does not influence much the value of the clusters polarizabilities. This comparison shows 
that the jellium model turns out to be quite a reasonable approximation. 

Figure O shows that the disagreement between theoretical and experimental points 
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is not always less than the experimental error bars. Such a disagreement might indicate 
that for certain N there have been experimentally detected cluster isomers other than 
those calculated in our work. For example, the calculated value = 659a. u. lies 

beyond the experimental error bars, while a N 3 £ g = 706. 876a. u. is within the range of the 
experimental error. 

Note that the polarizabilities of clusters Na 8 , Na w and Na 2 o, possessing the T d 
point symmetry group, surpass a little the corresponding experimental values, being 
quite close to them. For the Nag and Naio clusters, the disagreement of the theoretical 
and experimental values is within the range of the experimental error. The similar 
situation occurs for the Nan cluster, characterized by the C 2v point symmetry group. 
This cluster likely belongs to the cluster chain leading to the formation of the tetrahedron 
Na 2 Q cluster from the tetrahedron Na$ one (see our discussion in section |3~l| ). Such a 
situation allows us to assume that the polarizabilities of other clusters of this chain, 
which we have not analized in this paper, because they are energetically not favorable, 
will be also quite close to the experiment. 

3.6. Normal vibration modes 

Using the B3LYP method, we have calculated the normal vibration frequencies for the 
optimized neutral sodium clusters. The results of this calculation are shown in figure 
|15| . In this figure, we indicate the point symmetry group for those clusters for which 
more than one cluster isomer has been considered (see figure |l|) . Numerous frequencies 
shown in figure [15] are degenerate or nearly degenerate. This explains why the total 
number of frequencies for most of clusters is less than the number of vibrational degrees 
of freedom available in the system. In the more symmetric clusters, like Na-?, Na$, Na\o 
or Na 20 , the rate of generacy of the normal vibration modes is higher. 



Structure and properties of small sodium clusters. 



34 




Figure 15. Normal vibration frequencies calculated by the B3LYP method for the 
neutral sodium clusters with N < 20. For each cluster we mark the breathing mode 
in the spectrum by dotted line and the surface quadrupole vibration modes by dashed 
lines. The number near some of the lines indicate the degeneracy of the corresponding 
mode. Note that we make this only for quadrupole surface vibration modes. 



Knowledge of normal vibration modes and their frequencies is important for 
physical understanding and quantitative description of the relaxation of electron 
plasmon excitations in metal clusters [42]. One can visualize normal vibration modes, 
showing the directions and amplitudes of the atoms displacements by corresponding 
vectors. Since it is difficult to show all such pictures in this paper due to their large 
number. We focus instead only on the two types of modes breathing and quadrupole 
surface vibration modes. Namely these modes have been considered in |42| within the 
dynamical jellium model |^T[ for the treatment of the electron-phonon coupling in the 
spherical metal clusters Na2o, Na^o arid Na$2- 

In this paper, we discuss the appearance of these specific vibration modes in a 
cluster system and compare their frequencies with the predictions made in [|2] on the 
basis of the jellium model. For this purpose, we have analysed all calculated vibration 
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Figure 16. Surface and volume vibration modes for the selected neutral sodium 
clusters. Number near each cluster image indicates the frequency of the corresponding 
normal vibration mode. The values are given in cm -1 . 



modes and identified the breathing and three quadrupole vibrations for each cluster. In 
figure [16|, we present images of the breathing and quadrupole vibration modes for some 
clusters to illustrate the way, how the identification of the modes has been performed. 
This figure shows that the identification made is definite enough. 

The results of this analysis are shown in figure [15|, where for each cluster we mark 
the breathing mode in the spectrum by dotted line and the surface quadrupole vibration 
modes by dashed lines. The number near some of the lines indicate the degeneracy of 
the corresponding modes. Note that we make this only for quadrupole surface vibration 
modes. The degeneracy rate and the number of quadrupole surface vibration modes 
can be easily understood with the help of the cluster images shown in figure |l|. This 
figure shows that the prototype of the breathing mode exists already in the Na% and 
iVa 4 clusters. For the Na^ cluster, one can identify the quadrupole surface vibration 
mode, although it is meaningful to discuss surface vibrations only for the Na§ cluster 
and larger. Figure shows the frequencies of the breathing and surface vibration 
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modes decrease systematically with increasing cluster size, although this decrease has 
numerous irregularities, particularly for the clusters with N < 8. The frequency of 
the breathing mode decreases faster with the growth of N than the frequency of the 
quadrupole surface vibration mode. 

Let us compare the calculated frequencies of the breathing and surface vibration 
modes with the predictions of the jellium model. In jllj, it was shown that the 
breathing vibration mode frequencies calculated for the spherical Na 2Q , Na^ and Na 92 
respectively within the framework of the dynamical jellium model are quite close to the 
values derived from the phonon dispersion law for metals |57| 

2 _ 6VpK 

~ M Na (9 + k 2 v 2 F rl) ' {Zb) 

where M^ a = 4.2 • 10 4 is the mass of sodium atom, Vp = (97r/4) 1//3 /r is the velocity 
of cluster electrons on the Fermi surface, ro is the Wigner-Seitz radius. In the long 
wave limit, equation (|25|) reduces to the Bohm-Staver formula for the velocity of sound, 
dVt/dk = vp/ \/3MN a m 3 • 10 5 cm/s. This number is quite close to the real value of the 
velocity of sound in the bulk sodium: 3.2 • 10 5 cm/s. 

Using the dispersion low (p5|) , we estimate the breathing mode frequencies for 
the magic Na$ and Na 2 Q clusters. The results of this calculation are as follows 
^Na 8 — 104.09cm -1 , Qnci2 = 80.49cm -1 . In this calculation we have used r = 4. 

The frequency values obtained from (^) are close to those presented in figure pl| 
^Nas — 127.15cm -1 , f2jva 20 = 78.11cm -1 . The agreement of the frequencies is rather 
good for the Na 2 o cluster case. For Na$, the agreement is reasonable, but not as good 
as for Na 20 . Some disagreement arises due to the fact that the Wigner-Seitz radius for 
the Nag cluster is about 10% smaller than its bulk value. Indeed, substituting ro = 3.6 
in (^) one derives flNa 8 = 127.10cm -1 , which is in the nearly perfect agreement with 
the ab initio result. The decrease of the Wigner-Seitz radius can be easily understood 
from the analysis of the cluster geometry shown in figure [I]. 

Now let us compare the quadrupole surface vibration mode frequencies calculated 
in our paper (see figure [15]) with those following from the dynamical jellium model. 
According to ||42|| , the quadrupole surface vibration frequencies, Q 2 , for the spherical 
Na 2 o, Na^Q and Nag 2 clusters are equal to 56.48cm -1 , 48.41cm -1 and 32.28cm -1 , 
respectively. The value of the quadrupole surface vibration frequency for the Na 2 o 
cluster calculated in the present work is equal to 63.15cm -1 , which is rather close to the 
value predicted in [Q. 

The values of the quadrupole surface vibration frequencies calculated for Na 20 , Na i0 
and iVag2 show relatively slow decrease with the growth cluster size. Extrapolating these 
values towards smaller cluster sizes, we derive frequency values, which are consistent 
with those shown in figure pl| This comparison demonstrates that the jellium model 
calculation of the surface vibration frequencies is in a reasonable agreement with the 
more accurate ab initio many-body theory. 

The comparison of the jellium model results with those derived by the more accurate 
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ah initio many-body theory is important, because it forms theoretical background for 
the jellium model calculations in larger cluster systems, for which ah initio methods are 
hardly possible. The comparison with the jellium model, which we performed in this 
paper, can be extended towards larger cluster sizes and other collective modes of ions 
motion. 

4. Conclusion 

In this paper we have calculated the optimized structure and various characteristics of 
sodium clusters consisting of up to 20 atoms. We have used three different methods: 
B3LYP, MP4 and HF. It was demonstrated that the first two methods due to 
accounting for many-electron correlations provide much better agreement with the 
available experimental data and theoretical results based of the configuration interaction 
method as compared to that for the Hartree-Fock approximation. This was checked for 
various cluster characteristics: cluster geometries, binding energies per atom and the 
ionization potentials. 

We have also calculated and analyzed the dependence of the ionic component 
and total quadrupole moments of sodium clusters as a function of their size. It was 
demonstrated that the cluster shapes characterized by the quadrupole moments are in 
a reasonable agreement with the predictions of the jellium model and the results of the 
experimental observations. 

We have determined the normal vibration modes and their frequencies for a number 
of clusters and demonstrated their qualitative agreement with the predictions based on 
the jellium model. 

The results of this work can be extended in various directions. One can use 
the similar methods to study structure and properties of various types of clusters. 
It is interesting to extend calculations towards larger cluster sizes and perform more 
comparison with the results following from the jellium model and other simplified 
theories, based either on pseudopotentials or effective interatomic potentials. A lot of 
novel problems arise, when considering collisions and electron excitations in the clusters 
with the optimized geometries. These and many more other problems on atomic cluster 
physics can be tackled with the use of methods considered in our work. 
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Appendix A. Tables 

In Appendix, we present tables of the essential cluster characteristics. The binding 
energies per atom for neutral and singly charged clusters are compiled in tables [XT] 
and [A2| . The principal values of the quadrupole moment tensor for neutral and singly 
charged clusters are presented in tables 1X3] and 
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N 


Symmetry 


Eat (a.u.) 






HF/6-31 1 H/TH 

1 1.1 / U Oil \J \ LI 1 


MP4/6-31 1 d((] n) 

1V11 H: / U Oil VI 1 U ; |J 1 


R3T.YP/6-31 1 fifVn 


Rpf 1fil 


1 




-161.8459 


-161.8459 


-162.2866 




2 


Dooh 


393 RQ1 1 


393 71 4Q 


394 ^QQQ 


ooo oi 7/2 
-OZO.Ol i U 


3 


J^ooh 

Co 
b.C-2v 
^ Ah 


48^ ^/LD^ 

-485.5403 
-48^ ^38^ 

t: (J O . O (J <J 

-485 5282 


48^ ^626 

-485.5653 
-48 ^ ^6^6 

t: (J tj . <J \J iJ \J 

-485 5626 


486 8Q63 

-486.8960 
-486 8939 

t: (J VJ . (J C/ O iy 

-486 8889 


-484.9729 


4 


- L -'2/i 

Ddh 


-647 3871 
-647 3897 


-647 4433 


-649 2076 
-649 1965 


-646 6494 


5 




-809.2518 


-809.3008 


-811.5164 


-808.3174 


6 


D 3h 


-971.0915 
-971.0998 


-971.1880 
-971.1872 


-973.8324 
-973.8344 


-969.9899 
-989.9884 


7 


D 5h 


-1 1 39 Q469 


-1 1 33 0634 


-1 1 36 1 430 

11<JU. 1t:OU 


-1 1 31 661 
-iioi.uuiu 


8 

o 


1 d 


-1294.8015 


-1294.9410 


-1298.4606 


-1293.3395 


9 


Civ 


1 4^6 6466 




1 46D 7^Q7 




1 


<^2 
Da,, 

Ca 

- 1 a 






1 693 07^8 
-1 693 0734 
-1623 0554 
-1623 0530 




11 


Ci 

w 1 






-1785 3737 
-1785 3726 




12 








-1947 6917 




13 


Ci 

w 1 






-2110.0045 




14 




- 


- 


-2272.3092 


- 


15 


c s 






-2434.6188 




16 


c. 






-2596.9370 




17 


c s 






-2759.2537 




18 








-2921.5704 




19 


A* 






-3083.8730 




20 








-3246.2015 
-3246.1981 





Table Al. In this table we present the total energies of the optimized neutral sodium 
clusters. Numbers of atoms in clusters are given in the first column. In the second 
column, the point symmetry groups of clusters are shown. In the next three columns, 
the cluster total energies derived by the HF, M P4 and B3LYP methods are compiled. 
For the sake of comparison, the total energies computed by the CI method in [jl6| are 
presented in the sixth column. 
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N 


Symmetry 


(a.u.) 






HF/6-311G(d) 
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-970 9749 
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7 




-1132.8278 
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-1135.9994 


8 


C2v 


-1 294 6866 
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-1 298 3082 


Q 


n„. 
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x w 
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c„ 
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19 
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20 
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21 


o h 






-3408.3434 



Table A2. In this table we present the total energies of the optimized singly charged 
sodium clusters. Numbers of atoms in clusters are given in the first column. In the 
second column, the point symmetry groups of clusters are shown. In the next three 
columns, the cluster total energies derived by the HF, MP4 and B3LYP methods are 
compiled. 
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Table A3. In this table we present the principal values of the quadrupole moment 
tensor calculated for neutral sodium clusters. The first column shows numbers of atoms 
in clusters. The second column gives their point symmetry groups. In the last three 
columns, the principal values Q xx , Q yy and Q zz are given. They have been computed 
by the B3LYP method. 
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Table A4. In this table we present the principal values of the quadrupole moment 
tensor calculated for singly-charged sodium clusters. The first column shows numbers 
of atoms in clusters. The second column gives their point symmetry groups. In the 
last three columns, the principal values Q xx , Q yy and Q zz are given. They have been 
computed by the B3LYP method. 



